######### This script will create realized distributions from the potential distributions created by MAXENT\

#### Establish base directory

base.dir = '/home1/99/jc152199/MAXENT/'
setwd(base.dir)

#### Load library

library('SDMTools')

#### Read in species realized distribution matrix

rdm = read.csv(paste(base.dir,'rdm.csv',sep=''))

#### Read in subregion ASCII and subregion codes

srs = read.asc('subregions.asc')



srs.codes = read.csv('srscodes.csv', header=T)

### Some rows in srs.codes are repeated, grab unique rows

srs.codes = unique(srs.codes)

##### Need to loop through ASCIIs for multiple species
##### First identify ASCII directory

ascii.dir = '/home1/99/jc152199/MAXENT/output/'

##### Identify species

spp = dir(ascii.dir)

#### Remove directories which aren't skinks

spp = spp[c(1,13:15,17:19)]

###### Now loop through and identify individual ASCII for each species
##### Don't forget, this needs to be done for both the microCLIM, expCLIM, and BIOCLIM distributions

for (s in spp)

	{
	
	#### Species loop
	
	for (m in c('microCLIM','BIOCLIM'))
	
		{
		
		#### Model loop
		
		#### Read in an ASCII
		
		p.asc = read.asc(paste(ascii.dir,s,'/',m,'/output/',s,'_',m,'.asc',sep=''))
		
		#### Create a blank ASCII
		
		base.asc = p.asc
		
		base.asc[which(is.na(base.asc)==F)]=0
		
		#### Read in the MaxEnt results for this species-model run and apply the threshold to remove very low suitability cells
		
		mr = read.csv(paste(ascii.dir,s,'/',m,'/output/maxentResults.csv',sep=''), header=T)
		
		#### Select the threshold (in this case the 'Balance Training Omission and Predicted Area Threshold') which is the most conservative threshold
		
		mrt = mr[,66]
		
		#### Apply the threshold
		
		p.asc[which(p.asc<mrt)]=0
		
		######### Identify the species line from the rdm
		
		s.rdm = rdm[which(rdm$spp_id==s),c(1:48)]
		
		#### Create a blank column in srs.codes
		
		srs.codes$occur = NA
		
		#### Match s.rdm with srs codes
		
		for (sr in unique(srs.codes$Subregion))
		
			{
			
			srs.codes$occur[which(srs.codes$Subregion==sr)]=s.rdm[1,which(names(s.rdm)==sr)]
			
			}
			
		#### Now get unique positions from p.asc

		s.pos = as.data.frame(which(is.finite(p.asc), arr.ind = T))

		# Add lat and long to this data frame of row/column positions

		s.pos$lat = getXYcoords(p.asc)$x[s.pos$row]
		s.pos$long = getXYcoords(p.asc)$y[s.pos$col]
		
		#### Extract data from the subregion ASCII
		
		s.pos$sr = extract.data(cbind(s.pos$lat,s.pos$long),srs)
		
		#### Remove rows from s.pos where srs.codes$occur equals zero
		
		pos2write = NULL
		
		#### Loop through subregions which have an occur value of 1 in srs.codes
		
		for (sr in srs.codes$Region_id[which(srs.codes$occur==1)])
			
			{
				
			#### Identify positions in subregions that are occupied
				
			t.pos = s.pos[which(s.pos$sr==sr),]
				
			#### Bind these positions into a dataframe
				
			pos2write=rbind(pos2write,t.pos)
				
			}
	
		### Close loop
			
		#### Extract data from p.asc using positions from pos2write
		
		pos2write$data2write = extract.data(cbind(pos2write$lat,pos2write$long),p.asc)
		
		#### Write data from pos2write$data2write back onto base.asc
		
		base.asc[cbind(pos2write$row,pos2write$col)]=pos2write$data2write
		
		#### Write out base.asc
		
		write.asc.gz(base.asc,file=paste(ascii.dir,s,'/',m,'/output/',s,'_',m,'_Realized',sep=''))

		}
		
	cat('\n',s,' - Completed\n',sep='')	
		
	}
	
#### Done

		
		



